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COMPUTATIONAL  APPROACHES  TO  THE  STUDY 
OF  REACTIVE  FLOW  MIXING 

I.  Introduction 

Our  understanding  and  eventual  ability  to  predict  the  complicated  interactions  occurring  in 
turbulent  reactive  flow  problems  is  imperative  for  many  combustion  modelling  applications. 
The  presence  of  turbulence  alters  mixing  and  reaction  times  and  heat  and  mass  transfer  rates 
which  in  turn  modify  the  local  and  global  dynamic  properties  of  the  system.  The  goal  of  this 
paper  is  to  first  describe  the  fundamental  physical  processes  we  must  model  in  studying  reactive 
flow  turbulence  and  then  to  show  how  detailed  modelling  of  these  phenomena  may  be  used  to 
help  construct  practical  phenomenological  turbulence  models  for  reactive  flows. 

If  the  full  set  of  conservation  equations  for  mass,  momentum,  and  energy  in  multispecies 
flow  could  be  solved  exactly  for  a  large  enough  range  of  time  and  space  scales,  turbulence  and 
turbulent  mixing  phenomena  would  be  contained  explicitly  in  the  solutions.  Unfortunately 
even  the  fastest  and  largest  computers  are  neither  fast  nor  large  enough  to  solve  a  complete 
problem  from  first  principles  even  though  the  fundamental  mathematical  model  is  adequate. 
The  basic  problem  in  doing  the  computation  arises  from  the  vast  disparity  in  time  and  space 
scales  among  the  fundamental  physical  processes.  For  example,  consider  a  system  in  which  the 
mean  free  path  X,  the  characteristic  scale  length  of  diffusive  mixing  in  reacting  shocks,  is  six 
orders  of  magnitude  smaller  than  the  macroscopic  convection  length,  L.  A  three-dimensional 
model  resolving  both  scales  would  need  ~106  computational  cells  in  each  spatial  dimension,  or 
a  total  of  — 1018  grid  points.  Perhaps  it  might  be  possible  to  use  106  grid  points,  but  certainly 
not  twelve  orders  of  magnitude  more.  In  trying  to  represent  all  of  the  macroscopic,  turbulent 
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mixing,  and  characteristic  molecular  dissipation  scales  in  a  single  calculation,  we  are  confronted 
with  perhaps  the  most  perplexing  multiple  time-scale  and  multiple  space-scale  problem  in  the 
detailed  modelling  of  reactive  flows.  We  must  develop  some  practical  algorithm  to  encompass 
these  widely  disparate  scales  in  the  same  calculation. 

In  modelling  turbulence  phenomena,  it  has  been  the  rule  to  try  to  develop  separate  equa¬ 
tions  which  describe  the  subgrid  scale  mixing  and  reactions.  These  are  then  coupled  to  the 
macroscopic  coarsely  resolved  fluid  dynamics  calculations  we  know  how  to  perform.  Here  the 
proper  separation  and  recoupling  procedures  for  subgrid  turbulent  mixing  and  macroscopic  flow 
are  not  nearly  as  obvious  as  the  analogous  problem  of  coupling  a  detailed  chemical  kinetics 
reaction  scheme  to  a  laminar  hydrodynamics  calculation  [1,2].  It  is  not  even  clear  how  to 
represent  some  of  the  important  transient  mixing  phenomena  qualitatively  let  alone  arrive  at 
quantitative  estimates  of  their  effect. 

In  Section  II  of  this  paper  we  discuss  the  fluid  dynamic  origins  of  turbulent  mixing  in 
order  to  obtain  an  intuitive  picture  of  the  range  of  phenomena  which  must  be  represented. 
This  information  is  then  used  in  Section  III  to  develop  a  list  of  resulting  requirements  which  we 
would  like  an  ideal  turbulent  mixing  submodel  to  have. 

First-principle,  "ab  initio'  calculations  are  discussed  in  Section  IV.  These  detailed  calcula¬ 
tions  provide  information  about  the  dynamic  characteristics  and  spectral  evolution  at  length 
scales  smaller  than  can  be  resolved  with  a  macroscopic  fluid  dynamics  model.  The  calculations 
are  often  extremely  expensive,  but  must  be  performed  to  provide  detailed  understanding  of 
microscale  processes.  Phenomenological  turbulence  models  cannot  be  expected  to  tell  us  any¬ 
thing  fundamental  about  turbulence  or  turbulent  mixing.  The  basic  physical  processes  must  be 
built  into  the  phenomenology,  they  cannot  be  derived  from  it.  Thus  a  series  of  idealized  ab  ini- 
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tio  calculations  may  be  used,  for  example,  to  delineate  the  hydrodynamic  channels  for  turbulent 
energy  reapportionment  and  provide  a  way  to  analyze  those  properties  which  must  be  built  into 
the  phenomenological  models. 

Section  V  describes  four  types  of  phenomenological  turbulence  models.  The  phenomeno¬ 
logical  modelling  problem  involves  the  representation,  not  the  resolution,  of  all  of  the  myriads 
of  individual  fluid  dynamic  degrees  of  freedom  within  each  macroscopic  cell.  We  seek  the  most 
efficient  representation  for  these  subgrid  degrees  of  freedom.  If  the  basic  processes  controlling 
a  particular  measurable  transient  or  turbulent  transport  effect  are  not  represented  in  a 
phenomenological  model,  certainly  the  model  cannot  predict  the  occurrence  of  the  effect. 
Section  V  is  concluded  with  a  description  of  a  dynamic  local  spectral  model,  which  would  allow 
us  to  incorporate  the  multiscale  and  transient  processes  appearing  specifically  in  reactive  flow 
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II.  The  Origins  of  Turbulence 

In  this  section  we  discuss  three  aspects  of  reactive  flow  mixing  which  when  considered 
together  give  a  consistent  picture  of  turbulent  flows.  By  thinking  in  terms  of  these  three 
different  aspects  and  their  interactions,  we  gain  an  inituitive  understanding  of  reactive  tur¬ 
bulence  and  the  ways  in  which  it  differs  from  the  classical  non-reactive  cases.  Figure  1  illus¬ 
trates  these  aspects  schematically.  The  first  aspect,  labelled  hydrodynamic  channels,  is  based  on 
the  idea  that  turbulence  begins  in  a  laminar  macroscopic  flow  as  the  onset  of  a  fluid  instability. 
The  energy  which  drives  this  macroscopic  instability  then  cascades  through  various  hydro- 
dynamic  channels  and  is  spread  convectively  over  a  broad  spectral  range.  By  concentrating  on 
the  fluid  instabilities  and  their  evolution  as  a  response  to  changing  fluid  conditions,  turbulent 
mixing  can  be  viewed  as  a  cascade  of  systems  of  distinct  instabilities  [3].  The  second  aspect 
focuses  on  the  idea  that  turbulent  mixing  is  dependent  on  the  contortion  of  reactive  surfaces 
which  originally  separate  Teactive  species.  We  must  first  follow  and  then  predict  the  behavior  of 
these  surfaces  as  they  move  and  stretch  with  the  fluid.  Steady  increase  of  the  reactive  surface 
area  enhances  molecular  mixing  and  speeds  reactions.  The  detailed  chemical  kinetics  and 
chemical  energy  released  by  these  reactions  is  the  third  aspect  of  the  reactive  flow  mixing  prob¬ 
lem.  Energy  release  in  gaseous  flows  promotes  expansion  and,  if  the  flow  is  sufficiently  exoth¬ 
ermic,  can  induce  buoyancy.  Thus  the  kinetics,  as  shown  in  Fig.  1,  feeds  back  into  the  hydro- 
dynamic  channels  and  the  cycle  is  closed. 

An  extensive  body  of  knowledge  and  valuable  expertise  has  been  developed  in  using 
models  of  turbulence  and  turbulent  mixing  in  constant  density  non-reacting  fluids  [4, 5,6,7]. 
However,  in  combustion  systems  there  are  interactions  which  modify  or  may  even  obviate 
these  models.  One  such  physical  mechanism  is  the  localized  change  in  the  rates  of  heat  and 
mass  transport  due  to  the  temperature  dependence  of  transport  coefficients.  Another  mechan- 
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ism  arises  from  the  strong  density  gradients  occurring  locally  in  the  combustion.  Non-linear 
phenomena  such  as  these  determine  the  turbulent  scale  lengths  and  eventual  chemical  reaction 
effects.  The  various  scale  lengths  covered  by  diffusive  mixing,  convection,  and  fluid  instabili¬ 
ties  are  bounded  by  a  dissipative  mean  free  path  scale  length  on  one  hand  and  the  macroscopic 
system  scale  on  the  other.  This  is  shown  schematically  in  Figure  2.  Since  energy  release  due 
to  chemical  reactions  can  only  occur  where  the  fuel  and  oxidizer  mix  molecularly,  turbulent 
combustion  is  driven  by  an  energy  source  which  fluctuates  on  the  scale  of  the  turbulent  mixing 
lengths.  This  situation  differs  markedly  from  the  classical  turbulent  flows. 

IIA.  Dynamic  Fluid  Instabilities  and  Hydrodynamic  Channels:  The  First  Aspect 

The  phenomena  described  above  are  examples  of  possible  hydrodynamic  channels  for  tur¬ 
bulent  energy  reapportionment.  Most  of  these  channels  involve  transient  effects  for  which 
detailed  balance  or  steady  state  arguments  are  inapplicable  [8],  They  involve  the  quick  onset  of 
a  relatively  short  wavelength  fluid  instability  which  results  from  a  slow  or  parametric  variation 
of  the  background  flow.  Thus  the  configuration  changes  from  one  of  local  hydrodynamic  stabil¬ 
ity  to  one  of  instability.  At  and  just  beyond  the  transition  point  of  marginal  stability,  the 
growth  rates  are  very  slow  and  the  unstable  wavelengths  relatively  long.  As  the  background 
flow  continues  to  evolve  into  the  unstable  regime,  maximum  growth  rates  increase  and  instabil¬ 
ity  spreads  over  a  broader  band  of  the  spectrum.  The  turbulence,  wich  is  the  non-linear  man¬ 
ifestation  of  these  instabilities,  arises  in  bursts  on  the  macroscopic  scale  and  is  an  intrinsically 
transient  phenomenon.  Thus  the  idea  of  a  steady  state  cascade  of  either  vorticity  or  turbulent 
energy  is  an  artificial  concept  because  so  many  of  the  channels  underlying  the  turbulent  process 
are  intermittent.  New  concepts  are  needed  to  describe  these  turbulent  bursts. 

Marginal  stability  methods  19,10,11]  can  play  a  very  useful  role  in  determining  the  tran¬ 
sient  interactions  between  macroscopic  flows  and  the  microscale  instabilities  which  they  induce. 
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In  splitter-plate  experiments,  for  example,  the  Kelvin-Helmholtz  instability  grows  initially  at 
wavelengths  characteristic  of  the  entrance  flow.  For  certain  experimental  designs  the  instability 
may  be  essentially  two-dimensional  [12,13],  A  vortex  street  begins  to  develop  and  vortices 
coalesce  with  rapidly  stretching  braids  of  interfacial  material  winding  up  in  the  cores.  Figure  10 
shows  a  numerical  calculation  of  such  an  instability.  When  the  local  Reynold's  number  is  high 
enough,  the  flow  in  the  vortices  becomes  Taylor  unstable  in  the  third  dimension  at  still  shorter 
wavelength.  Thus  a  one-time  situation  occurs  and  then  passes. 

In  splitter  plate  experiments  there  is  a  constant  source  of  energy  and  a  tendency  toward 
instability.  The  turbulence  grows  up  to  a  level  where  the  governing  flow  remains  near  marginal 
stability.  If  the  turbulence  decays  somewhat,  the  background  flow  becomes  more  unstable  until 
the  turbulence  builds  up  again.  If  the  turbulence  becomes  somewhat  too  strong,  the  back¬ 
ground  flow  is  broken  down  and  smeared  out  with  the  general  effect  of  reducing  the  growth  of 
instabilities.  In  the  smoother  flow  the  turbulence  soon  starts  to  abate.  Since  there  is  a  decay 
time  associated  with  the  turbulence  and  an  induction  time-lag  for  the  instabilities,  it  is  not 
surprising  that  intermittency  and  bursts  of  turbulence  are  the  norm  rather  than  the  exception. 

One  of  the  major  physical  conditions  which  distinguishes  reactive  flow  from  classical  tur¬ 
bulence  is  the  existence  of  density  gradients.  These  cause  the  baroclinic  generation  of  vorticity, 
£,  through  the  term 

d£  Vp  x  VP 

dt  ~  p  *  +"‘* 

where  p  and  P  are  the  multidimensional  fluid  density  and  pressure  fields,  respectively.  Below 
we  consider  two  of  the  hydrodynamic  channels  activated  by  density  gradients.  These  channels 
are  rather  typical  of  the  many  multi-scale  interactive  effects  which  will  have  to  be  understood 
before  we  can  claim  to  understand  turbulence. 
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When  the  mass  density  p  is  constant,  as  is  assumed  in  much  of  the  classical  turbulence 
literature,  all  of  the  effects  potentially  depending  on  the  vorticity  source  term  are  absent  by 
construction.  In  the  classical  picture,  the  turbulent  vorticity  spectrum  is  driven  at  macroscopic 
scales,  the  intermediate  wavelengths  are  populated  by  cascade  and  vortex  stretching,  and  the 
short  wavelengths  in  the  spectrum  are  dissipated  viscously  (collisionally).  In  combustion,  the 
localized  release  of  heat  in  molecularly  mixed  fuel-oxidizer  pockets  causes  strong  transient 
expansion  of  the  combusting  gases.  The  resulting  low  density  region  has  scale  lengths  charac¬ 
teristic  of  the  combustion  process  and  interacts  with  pressure  gradients  to  generate  vorticity  in 
the  flow  on  these  same  characteristic  scales.  Vorticity  on  these  scales  is  efficient  at  turbulent 
mixing  and  can  be  expected  to  feed  back  on  itself.  The  combustion  process  is  in  turn  enhanced 
as  long  as  fuel  and  oxidizer  are  present.  Since  the  spectrum  seems  now  to  be  driven  at  short 
wavelengths  as  well  as  long,  plateaus  or  even  peaks  might  form  in  the  spectrum  altering  the 
usual  notions  of  cascade,  scaling,  and  hence  modelling. 

There  are  two  different  channels  for  turbulent  mixing  that  this  fluid  expansion  activates  in 
chemically  reactive  flows.  These  can  be  understood  qualitatively  in  terms  of  localized 
Rayleigh-Taylor  growth.  The  first  role,  an  active  one,  occurs  when  a  pocket  of  gas  is  actively 
expanding  due  to  the  heat  released  from  chemical  reactions.  This  case  is  shown  schematically 
in  Figure  3.  The  expansion  occurs  in  a  restricted  volume  because  the  fuel  and  oxidizer  enter 
the  reaction  region  separately.  A  time-varying  acceleration  accompanying  this  expansion  is  felt 
in  the  surrounding  fluid.  A  region  of  strong  density  gradients  undergoing  this  acceleration  will 
be  subject  to  the  Rayleigh-Taylor  instability  during  much  of  the  time  the  energy  is  being 
released.  A  perturbation,  whose  wavelength  is  comparable  to  the  distance  the  fluid  moves 
while  expanding,  will  e-fold  about  once  from  the  expansion. 
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The  second  role  of  density  gradients,  depicted  in  Figure  4,  is  passive ,  but  is  potentially 
more  important  than  the  active  role.  In  the  passive  role  the  expansion  influences  turbulent  mix¬ 
ing  by  providing  the  density  gradient  which  leads  to  vorticity  generation.  Consider  a  localized 
fluid  vortex  rotating  at  angular  frequency  tu.  The  acceleration  of  the  fluid  is  Ru> 2  at  radius  R. 
Smaller  low-density  pockets  in  this  vortex,  formed  by  chemical  reaction  at  an  earlier  time  are 
unstable  to  the  Rayleigh-Taylor  modes  where  the  local  centrifugal  fluid  acceleration  points 
opposite  to  the  density  gradient.  The  acceleration  is  provided  by  the  divergence-free  rotational 
flow  generated  earlier  or  driven  into  the  system  externally. 

Because  the  active  generation  of  vorticity  is  limited  in  time  to  the  expansion  phase  itself, 
it  is  reasonable  to  expect  that  the  amount  of  vorticity  generated  by  the  passive  interaction  of 
existing  density  gradients  with  large  scale  vortices  is  much  greater.  A  very  simple  analysis  indi¬ 
cates  that  there  are  several  growth  times  of  the  Rayleigh-Taylor  mode  during  a  single  rotation 
of  the  vortex  when  the  mode  wavelength  is  comparable  to  the  vortex  radius.  Both  the  active 
and  the  passive  role  of  expansion  are  not  described  by  standard  incompressible  fluid  turbulence 
theories  and  models. 

IIB.  Reactive  Interface  Dynamics:  The  Second  Aspect 

Turbulent  mixing  can  also  be  viewed  as  the  stretching  and  convolution  of  surfaces  which 
originally  separate  reactive  species.  These  reactive  interfaces  or  surfaces  move  with  the  fluid. 
As  the  fuel  and  oxidizer  diffuse  molecularly,  these  surfaces  continue  to  lie  normal  to  the 
strongest  species  density  gradients.  The  situation  is  presented  in  Figure  5.  In  the  upper  left- 
hand  panel  a  cube  of  fluid  1  cm  on  a  side  is  shown  at  time  t  —  0.  An  interface  between  react¬ 
ing  species  A  and  B  divides  the  cube  into  two  halves.  Subsequent  stretching  of  this  reactive 
surface  is  controlled  by  the  local  fluid  velocity. 
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Motion  of  the  fluid  elements  making  up  the  reactive  surface  breaks  down  into  a  com¬ 
ponent  normal  to  the  surface  which  can  change  the  integrated  area  and  components  parallel  to 
the  surface  which  cannot.  All  of  these  components  are  important  and  complement  each  other 
in  the  turbulent  mixing  process.  When  nearby  points  on  the  surface  separate  or  stretch,  the 
fluid  on  opposite  sides  of  the  surface  approaches  the  surface  to  keep  the  flow  roughly 
divergence-free.  Any  species  or  temperature  gradients  normal  to  the  initial  surface  are 
enhanced,  which  in  turn  increase  the  diffusive  interpenetration  of  the  reactants.  This  surface 
stretching  process  is  independent  of  convective  interpenetration  of  material  which  also 
enhances  mixing.  This  latter  process,  shown  on  the  right  of  Figure  6,  is  governed  by  deriva¬ 
tives  of  the  normal  velocity  parallel  to  the  surface.  By  increasing  the  actual  area  of  the  reactive 
surface,  the  bulk  reactivity  is  also  increased.  Stretching  and  interpenetration  usually  occur 
simultaneously  but  for  clarity  have  been  illustrated  separately  in  Figure  6. 

If  the  turbulent  spectrum  is  dominated  by  short  wavelengths,  an  originally  smooth  surface 
will  become  very  wrinkled  at  short  wavelengths  before  larger  scale  convolutions  have  had  a 
chance  to  grow.  This  is  illustrated  in  the  upper  right  hand  panel  of  Figure  7.  Because  the 
length  scale  of  the  turbulent  spectrum  is  close  to  molecular  mixing  lengths,  the  effective 
volume  in  which  enhanced  mixing  occurs  is  much  less  than  what  would  be  calculated  knowing 
the  stretched  surface  area.  Adjacent  folds  of  the  surface  which  approach  within  the  molecular 
diffusion  length  of  each  other  tend  to  merge  their  mixed  volumes,  as  shown  in  the  Figure  7. 
The  small  surface  corrugations  rapidly  get  smoothed  over  with  diffusive  "fluff"  so  the  effect  of 
molecular  diffusion  is  to  limit  the  growth  of  the  reactive  surface  area. 

In  the  opposite  case,  shown  in  the  lower  left  hand  panel  of  Figure  7,  the  dominant 
wavelengths  are  long  compared  io  the  molecular  mixing  length.  The  surface  is  characterized  by 
relatively  smooth,  long  wavelength  bulges  which  can  fill  the  entire  cell  before  any  of  the  mutual 
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interference  effects  described  in  the  previous  paragraph  can  occur.  In  fact,  the  area  of  an  ideal 
Lagrangian  surface  in  an  isotropic,  homogeneous  turbulent  velocity  field  should  probably 
increase  exponentially.  Molecular  mixing  now  occurs  along  convolutions  which  have  scales 
longer  than  the  diffusion  length.  Eventually,  however,  overlap  of  the  effective  mixed  volume 
occurs  as  the  surface  becomes  more  convoluted.  Thus,  the  molecularly  mixed  volume  is  again 
limited  even  though  the  ideal  reactive  surface  area  goes  to  infinity. 

IIC.  Detailed  Chemical  Kinetics:  The  Third  Aspect 

The  two  aspects  of  time-dependent  turbulent  mixing  in  reactive  flows  presented  above, 
the  hierarchy  of  instabilities  and  the  evolution  of  the  reactive  surface,  focus  on  the  nonlinear 
fluid  dynamic  interactions  and  the  reactant-mixing  effects  respectively.  Thus  these  two  aspects 
complement  each  other  and  a  clear  picture  from  each  will  be  required  to  gain  a  proper  perspec¬ 
tive  on  turbulent  mixing.  The  third  aspect  we  wish  to  discuss  closes  the  schematic  loop  of  Fig. 
1.  The  essentially  one-dimensional  profiles  of  the  fuel,  oxidizer,  and  other  reactants  perpendic¬ 
ular  to  the  reactive  surface  change  not  only  from  convection  and  molecular  diffusion  but  also 
from  chemical  reactions.  The  fluid  dynamic  expansion  and  heat  release  must  self-consistently 
couple  back  into  the  excitation  of  hydrodynamic  channels  to  complete  the  picture. 

A  major  concern  then  is  the  range  of  validity  of  fast  chemistry  approximations  made  to 
simplify  the  analysis.  In  some  cases  the  overall  reaction  rate  is  governed  by  the  diffusion  of  hot 
fuel  and  oxidizer  together  through  an  expanding  region  of  hot  products  and  reactants.  In  other 
cases  the  fuel  and  oxidizer  mix  convectively  and  then  molecularly  before  ignition  occurs. 
These  latter  situations  are  generally  ignited  by  a  rapidly  moving  flame  front  which  travels  paral¬ 
lel  to  the  reactive  interface  rather  than  perpendicular  to  it.  Here  finite  chemical  kinetics  clearly 
plays  a  crucial  role  since  it  determines  the  flame  speed. 
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Figure  1  emphasizes  the  feedback  and  interactions  between  the  study  of  hydrodynamic 
channels,  of  reactive  interface  dynamics,  and  of  detailed  chemical  kinetics.  These  three  distinct 
aspects  also  suggest  independent  yet  complementary  approaches  to  ab  initio  computational  tech¬ 
niques,  a  subject  considered  further  in  Section  IV. 
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111.  Sume  Properties  of  an  "Ideal"  Subgrid  Mixing  Model 

Accurate  yet  compact  phenomenological  turbulence  models  must  be  developed  to  model 
realistic  combustor  systems,  open  flames,  and  other  turbulent  reactive  flows  confidently  and 
efficiently.  These  computational  models  must  asymptotically  decouple  ihe  subgrid  turbulence 
and  microscopic  instability  mechanisms  from  calculations  of  the  macroscopic  flow.  In  this  sec¬ 
tion,  the  important  properties  which  we  would  like  to  see  incorporated  in  an  ideal  turbulence 
model  are  listed  and  described.  The  goal  is  to  use  these  desired  properties  to  guide  us  in  both 
extracting  information  from  ab  initio  calculations  (Section  IV)  and  constructing  phenomenolog¬ 
ical  turbulence  models  (Section  V). 

Ill  A.  Chemistry  —  Hydrodynamic  Coupling  and  Feedback 

Explicit  energy  feedback  mechanisms  from  the  mixing  and  reactions  to  the  turbulent  velo¬ 
city  field  and  the  macroscopic  flow  must  be  formulated.  The  "laminar"  macroscopic  flow  equa¬ 
tions  contain  phenomenological  terms  which  represent  averages  over  the  macroscopic  dynamics. 
Examples  of  these  terms  are  eddy  viscosity  and  diffusi vity  coefficients  and  average  chemical 
heat  release  terms  which  appear  as  sources  in  the  macroscopic  flow  equations  These  modified 
Navier-Stokes  equations  are  postulated  to  include  the  effects  of  turbulence.  Besides  providing 
these  phenomenological  terms,  the  turbulence  model  must  make  use  of  the  information  pro¬ 
vided  by  the  large  scale  flow  dynamics  in  order  to  determine  the  energy  which  drives  the  tur¬ 
bulence.  The  model  must  be  able  to  follow  reactive  interfaces  on  the  macroscopic  scale. 

IIIB.  Modelling  Onset  and  Other  Transient  Turbulence  Phenomena 

The  model  should  be  able  to  predict  the  onset  of  turbulence  in  what  was  initially  laminar 
flow  since  bursts  and  other  highly  transient  phenomena  seem  to  be  the  rule  in  reactive  flow  tur¬ 
bulence.  The  fundamental  gradients  in  density,  temperature,  and  velocity  fields  in  the  reacting 
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fluid  drive  the  macroscopic  fluid  dynamic  instabilities  which  initiate  turbulence.  Thus  these  gra¬ 
dients  from  the  macroscopic  calculation  are  bound  to  be  key  ingredients  in  determining  the 
energy  that  is  available  to  drive  the  turbulence.  Density  stratification  in  a  time-dependent  fluid 
dynamics  model  is  essential. 

IIIC.  Complicated  Reactions  and  Flow 

The  ideal  turbulence  model  should  allow  for  detailed  calculations  or  suitable  parameteriza- 
tions  of  chemical  kinetics,  of  buoyancy  effects,  and  of  the  other  hydrodynamic  channels  which 
the  physics  in  a  given  situation  might  require.  In  particular  it  must  be  possible  to  deal  with 
multiscale  effects  within  the  subgrid  model.  If  there  is  a  delay  as  velocity  cascades  to  the  short 
wavelength  end  of  the  spectrum,  the  model  must  be  capable  of  representing  this.  Otherwise 
bursts  and  intermittency  phenomena  cannot  be  calculated. 

IIID.  Lagrangian  Framework 

An  ideal  subgrid  model  should  be  constructed  on  a  Lagrangian  hydrodynamics  framework 
moving  with  the  macroscopic  flow.  This  requirement  reduces  purely  numerical  diffusion  to 
zero  so  that  realistic  turbulence  and  molecular  mixing  phenomena  will  not  be  masked  by  non¬ 
physical  numerical  smoothing.  This  requirement  also  removes  from  concern  the  possibility  of 
masking  purely  local  fluctuations  by  truncation  errors  occuring  when  macroscopic  convective 
derivatives  are  represented  numerically.  The  time-dependent  (hyperbolic)  Lagrangian  frame¬ 
work  should  also  generalize  to  three  dimensions  as  well  as  resolve  reactive  interfaces  dynami¬ 
cally. 

IIIE.  Scaling 

As  mentioned  in  Section  II,  breaking  the  calculation  into  macroscopic  scales  and  subgrid- 
scales  is  an  artifice  to  allow  us  to  model  turbulence.  The  important  physics  occurs  continuously 
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over  the  whole  spectrum  from  k  —  0  to  kdlss,  which  is  the  wave  number  corresponding  to  a 
mean  free  path  of  a  molecule.  Thus  the  macroscopic  and  subgrid  scale  spectra  of  any  physical 
quantity  must  couple  smoothly  at  kce„  the  cell  boundary  wave  number.  If  this  number  were  to 
be  changed,  as  might  happen  if  numerical  resolution  were  halved  or  doubled,  the  predictions  of 
the  turbulence  model  must  not  change. 

IIIF.  Efficiency 

Of  course,  the  model  must  be  efficient.  The  number  of  degrees  of  freedom  required  to 
specify  the  status  of  turbulence  in  each  separately  resolved  subgrid  region  has  to  be  kept  to  a 
minimum  for  the  model  to  be  generally  useable.  The  real  fluid  has  essentially  an  infinite 
number  of  degrees  of  freedom  to  represent  the  state  of  the  gas  in  each  small  element.  We 
would  like  to  be  able  to  do  the  job  with  a  minimal  number  of  degrees  of  freedom. 

How  to  choose  a  representation  to  have  these  ideal  properties  is  not  obvious.  Several 
alternative  approaches  are  mentioned  in  Section  V  of  this  paper.  It  is  obvious  that  considera¬ 
tion  of  the  dynamic  aspects  of  instability  cascade  and  reactive  surface  evolution  necessitate 
significant  departures  in  modelling  techniques  from  methods  found  adequate  when  fast  chemis¬ 
try  or  steady  state  Eulerian  models  of  classical  turbulence  are  assumed. 
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IV.  Ab  Initio  Calculations 

In  order  to  assess  the  interaction  between  energy  release  and  turbulence  quantitatively,  it 
is  possible  to  perform  idealized  fundamental  reactive  fluid  dynamic  calculations  taking  advan¬ 
tage  of  today’s  advanced  computer  systems.  Available  practical  algorithms  permit  us  to  solve 
the  full  set  of  conservation  equations  for  multi-dimensional  multi-species  flow  on  what  is  gen¬ 
erally  considered  the  hydrodynamic  subgrid  scale.  These  numerical  models  consider  buoyancy 
and  compressibility  effects  and  include  realistic  expressions  for  heat  and  mass  diffusion  and 
chemical  reactions.  Thus  the  full  spectrum  of  compressible  fluid  dynamic  effects  may  be  simu¬ 
lated  in  a  small  region  where  enough  resolution  is  available  to  resolve  the  turbulence  and  the 
dissipation  scales  as  well  as  the  necessary  macroscopic  convective  effects. 

The  ultimate  goal  of  such  "ab  initio"  calculations  is  to  advance  our  basic  understanding  of 
turbulence  and  to  provide  information  about  the  small  scale  mixing  and  energy  release  which 
can  then  be  used  to  construct  and  calibrate  phenomenological  turbulence  models.  In  constant 
density,  idealized,  incompressible  turbulence,  a  number  of  these  ab  initio  calculations  have 
been  successfully  performed  [14].  Such  calculations  could  also  be  made  for  compressible  reac¬ 
tive  turbulence,  but  the  computations  are  somewhat  more  difficult.  This  difficulty  manifests 
itself  in  a  number  of  ways.  First,  resolution  of  all  scales  is  not  generally  possible  so  only  very 
idealized  and  restricted  questions  may  be  asked.  Second,  an  ab  initio  calculation  is  performed 
only  for  one  specific  initial  condition.  Thus  statistical  information  about  the  "turbulence"  has  to 
be  obtained  in  transient  problems  by  performing  a  number  of  calculations  and  averaging. 
Finally,  the  cost  of  these  detailed,  first-principle  fluid  dynamic  calculations  is  usually  high. 

The  results  of  these  detailed  ab  initio  simulations  are  correspondingly  valuable.  Any 
correlation,  stress  term,  fluctuation,  or  overlap  integral  can  be  determined  from  the  computed 
flow  without  interfering  with  the  numerical  experiment.  Therefore  closure  assumptions  for  the 
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more  phenomenological  decoupled  turbulence  models  discussed  in  Section  V  may  be  tested 
directly.  This  calibration  procedure  is  particularly  valuable  because  the  exact  computational 
solution  of  many  important  large-scale  reactive  (low  problems  is  far  beyond  present  technology. 
It  will  be  necessary  to  use  intermediate  phenomenological,  i.e.,  lumped-parameter,  models  to 
incorporate  all  of  the  scales  for  practical  engineering  applications.  These  models  must  be  care¬ 
fully  calibrated  using  experimental  data  and  specific  ab  initio  calculations. 

The  ab  initio  calculations  themselves,  not  their  parameterizations,  are  the  only  computa¬ 
tions  that  can  teach  us  something  about  subgrid  scale  turbulence.  Realistic  chemical  models 
can  be  included  on  the  microscopic  scale  if  a  small  enough  system  is  being  treated.  Therefore 
each  fluid  element,  as  it  undergoes  hydrodynamic  motions  which  differ  slightly  from  its  neigh¬ 
bors,  experiences  a  correspondingly  different  chemical  history.  The  lumped-parameter  subgrid 
model  must  necessarily  represent  an  average  behavior  in  a  macroscopic  fluid  element.  These 
averages  can  be  evaluated  directly  in  idealized  situations  using  large  detailed  ab  initio  tur¬ 
bulence  simulations. 

The  probability  density  functions  (PDF’s)  currently  used  in  multi-moment  modelling  of 
turbulence  can  be  measured  from  detailed  simulations  as  a  function  of  time,  configuration,  and 
chemical  kinetics.  This  can  even  be  done  in  situations  where  the  energetic  feedback  of  the 
reactions  on  the  hydrodynamics  is  specifically  incorporated.  Thus  not  only  can  phenomenologi¬ 
cal  models  be  tested,  they  can  be  upgraded  as  improved  phenomenologies  are  derived  from  the 
analyses  of  the  ab  initio  calculations  themselves. 

Perhaps  most  important,  ab  initio  calculations  can  be  used  to  answer  fundamental  ques¬ 
tions  about  the  turbulent  mixing  process  itself.  In  Section  II,  the  active  and  passive  roles  of 
expansion  in  reactive  flow  turbulence  were  discussed  as  an  illustration  of  nonlinear  fluid 
dynamic  phenomena  which  enter  when  sharp  density  variations  and  energy-significant  reactions 
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characterize  the  flow.  These  phenomena  are  examples  of  many  identifiable  hydrodynamic  chan¬ 
nels  for  turbulent  energy  reapportionment.  Ab  initio  calculations  can  be  used  to  isolate  these 
individual  channels  and  study  their  characteristic  signature  in  the  flow.  Other  examples  of  such 
channels  would  be  mixing  from  growth  of  the  Kelvin-Helmholtz  instability  at  shearing 
interfaces  driven  by  the  effects  of  buoyancy,  Taylor  instability  arising  naturally  in  larger  vortices 
and  vortex  streets,  and  Rayleigh-Taylor  instability  mixing  where  the  pressure  gradients  are 
gravitationally  induced.  Each  of  these  channels  is  really  a  mechanism  for  taking  relatively 
organized  kinetic,  potential,  or  internal  energy  on  one  scale  and  reapportioning  it  to  other 
scales. 

Results  from  one  such  ab  initio  calculation  are  shown  in  Figures  8  and  9.  A  two- 
dimensional  version  of  the  FCT  reactive  flow  model  [2,15]  was  used  to  solve  the  conservation 
equations  for  mass,  momentum  and  energy  on  a  square,  doubly  periodic  domain  using  an  ideal 
gas  equation  of  state.  A  three  species  model  kinetics  scheme  is  included  to  represent  the 
exothermic  reaction  A  +  B  —  C  +  Heat.  Each  of  the  three  species  mass  densities  (pA,  pB, 
and  pc)  satisfies  its  own  continuity  equation  with  an  added  molecular  diffusion  term.  Viscosity 
is  neglected  (a  von  Neuman-like  treatment  of  shock  heating  is  not  needed  in  the  FCT  formal¬ 
ism).  Figures  8  and  9  show  different  computer  plots  of  the  results  for  two  fluids  of  different 
density,  which  interpenetrate  and  diffusively  mix  as  the  result  of  an  initial  velocity  field.  Using 
such  calculations  we  hope  to  measure  the  variation  of  mixing  rate  with  density  difference. 
From  this  calculation  we  see  that  "ab  initio"  models  can  only  be  used  in  very  idealized  tur¬ 
bulence  problems  although  these  limited  problems  can  be  solved. 

Another  two-dimensional  ab  initio  calculation  performed  on  the  Kelvin-Helmholtz  insta¬ 
bility  [16]  is  shown  in  the  four  panels  of  Figure  10.  Two  different  fluids  initially  in  relative 
shear  flow  abut  at  a  sharp  interface  at  x  -  0.  The  fluid  on  the  left  moving  up  and  the  fluid  on 
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the  right  moving  down.  Solid  wall  boundary  conditions  are  applied  to  the  left  and  right  boun¬ 
daries  at  \  —  +  and  x  “  —respectively.  The  vertical  boundaries  are  periodic  to  accommodate 
the  inflow  and  outflow  systematically.  Enough  cells  are  present  in  the  calculation  (200  x  360)  to 
allow  representation  of  relatively  short  scale  lengths  such  as  the  interface  between  species  A 
and  B  (shown  as  two  contours  at  the  levels  0.4,  and  0.6  in  the  function  paKpa  +  p«).  Inter¬ 
mediate  scale  lengths  such  as  the  vortex  size  which  develops  and  longer  scales  characteristic  of 
the  system  size  can  thus  be  easily  resolved.  The  interaction  of  these  scale  lengths  can  be 
displayed  and  studied  explicitly  in  properly  posed  and  analyzed  ab  initio  calculations. 

Two  additional  data  curves  superimposed  on  the  contour  plots  of  Figure  10  are  the  one¬ 
dimensional  vertical  averages  of  the  species  densities  at  each  horizontal  location  in  the  calcula¬ 
tion.  These  averages  appear  as  solid  bars  at  y  —  0  ant  y  —  in  the  upper  left  hand  panel,  our 
initial  configuration.  Species  A,  plotted  as  A’s,  has  a  vertically  averaged  value  of  pA/p  —  1  on 
the  left  of  the  initial  shear  interface  and  zero  on  the  right.  Species  B  has  just  the  reverse  initial 
profile.  As  the  fluid  instability  progresses,  the  vertically  averaged  density  shows  a  gradual  tansi- 
tion  crossing  the  shear  layer,  even  though  the  spatially  varying  and  stretching  reactive  surface 
remains  sharp.  In  this  way,  deterministic,  detailed,  convective  effects  in  the  ab  initio  calcula¬ 
tions  can  be  compared  one-for-one  with  fuzzy,  statistical,  averaged,  diffusive-like  effects  in  a 
phenomenological  model. 

Also  from  Figure  10  we  can  see  that  the  average  density  displays  a  rather  discontinuous 
drop  at  the  edge  of  the  turbulent  layer  reminiscent  more  of  erosion  at  the  edge  of  a  cliff  than  of 
gentle  interdiffusive  penetration  of  two  species.  The  effective  one-dimensional  diffusion 
coefficient  has  to  be  highly  variable  in  space  and  time  to  model  the  intrinsically  two-dimensional 
mixing  which  occurs.  Such  variability  may  be  provided  by  making  the  turbulent  diffusivity  a 
nonlinear  function  of  the  "turbulence"  itself.  We  can  also  see,  particularly  in  the  lower  right 
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hand  panel,  that  even  though  both  <p4/p>  and  <pB/p>  are  about  1/2  at  x  “  0,  very  little 
of  the  fluid  is  actually  mixed.  Rather,  strips  of  one  or  the  other  pure  fluid  are  interleaved. 
Reactions  between  A  and  B  can  only  occur  where  they  are  molecularly  mixed.  Thus  reactions 
would  still  be  occurring  as  a  surface  phenomena  even  though  <pA>  <pB>  in  an 
apparently  extensive  volume. 

This  last  example  makes  clear  one  of  the  great  operational  advantages  of  using  detailed  ab 
initio  calculations  in  the  construction  of  phenomenological  turbulence  models.  A  crudely 
resolved  macroscopic  flow  calculation  with  a  subgrid  model  for  turbulence  could  be  compared  to 
an  ab  initio  calculation  of  the  same  exact  problem.  By  varying  only  the  resolution  and  keeping 
the  initial  conditions,  fluid  dynamics  algorithm,  and  boundary  conditions  identical,  the  transi¬ 
tion  from  an  ab  initio  to  a  phenomenological  representation  can  be  studied  in  detail. 

In  a  complete  model  of  reactive  turbulence,  all  space  and  time  scales  between  the  macros¬ 
copic  convection  and  the  molecular  dissipation  scales  are  equally  important  to  the  fluid  and 
must  enter  into  the  description  of  turbulence.  There  is  no  intermediate  scale  which  can  be 
ignored  safely  in  making  the  asymptotic  decoupling  between  macroscopic  and  statistically  indis¬ 
tinguishable  microscopic  scales.  However,  the  transition  from  a  discrete  representation  on  a 
macroscopic  mesh  to  a  phenomenological  subgrid  representation  is  necessarily  discontinuous 
and  therefore  likely  to  be  a  major  source  of  difficulty. 

In  the  ab  initio  calculations,  the  characteristics  of  and  interactions  between  the  various 
hydrodynamic  channels  for  turbulent  reapportionment  comprise  a  local  description  of  the  tur¬ 
bulent  mixing  problem.  This  local  approach  risks  missing  global  interactions  in  and  constraints 
on  the  flow.  These  can  only  be  treated  by  using  a  complementary  phenomenological  model 
coupled  with  the  large  scale  macroscopic  flow.  Various  options  for  phenomenological  modelling 
are  discussed  in  the  next  section. 
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V.  Options  for  Phenomenological  Turbulence  Modelling 

Most  of  the  current  “phenomenological*  turbulence  modelling  techniques  are  based  on  the 
idea  that  local  fluctuations  in  the  macroscopic  mean  flow  decorrelate  rapidly  in  space.  Thus  the 
local  turbulence  details  would  depend  only  weakly  on  turbulent  fluctuations  a  small  but  macros¬ 
copic  distance  away.  Within  this  context,  a  local  turbulent  energy  density,  a  characteristic  mix¬ 
ing  length,  and  concentration  fluctuations  may  be  defined.  Conversely,  fluctuating  coherent 
structures  on  the  macroscale  must  be  represented  in  the  fluid  dynamics;  a  local  turbulence 
model  cannot  be  expected  to  deal  with  non-local  fluid  dynamics. 

Below  we  describe  four  kinds  of  phenomenological  models: 

1.  Turbulent  Energy  and  Scale  Models, 

2.  Moment  Equation  Methods, 

3.  Physically  Motivated  Phenomenological  Models 

4.  Localized  Spectral  Dynamics  Model. 

These  approaches  are  discussed  in  terms  of  the  six  criteria  listed  in  Section  III.  The  four  gen¬ 
eric  model  types  described  below  are  certainly  not  all  inclusive,  but  they  are  representative  of 
the  major  approaches  pursued  today  in  reactive  turbulence. 

VA.  Turbulent  Energy  and  Scale  Models 

Based  on  the  idea  that  turbulence  can  be  described  by  local  fluctuations,  a  number  of 
authors  have  written  dowi  evolution  equations  to  propagate  the  turbulence  energy,  scale 
lengths,  species  overlap  integrals,  etc.  The  equations  for  these  subgrid  quantities  or  combina¬ 
tions  of  them  take  the  same  form  as  those  of  the  mean  fluid  quantities.  The  characteristic 
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source  and  sink  terms  are  functions  of  both  mean  flow  variables  and  local  turbulence  parame¬ 
ters.  There  is  a  long  history  of  work  in  this  area,  starting  with  Kolmogorov  (17)  in  1942  and 
continuing  to  date  with  the  three-equation  model  of  Spalding  [18]. 

These  one-,  two-,  and  three-equation  approaches  have  been  used  most  successfully  for 
certain  classical  non-reactive  flow  problems  which  primarily  involve  shear  flows  with  high  Rey¬ 
nolds  numbers  away  from  boundaries  and  turbulent  boundary  layers.  The  results  are  often 
acceptably  accurate  in  describing  problems  in  which  the  turbulence  has  already  developed. 
Various  free  parameters  must  be  optimized,  or  calibrated  against  experiments.  These  methods 
are  questionable  for  the  treatment  of,  for  example,  reactive  flows,  high  mach  number  flows, 
low  Reynolds  number  flows,  flows  where  gravity  is  important,  or  intermittant  turbulent  bursts. 
The  presence  of  large  coherent  structures  in  real  mixing  situations  casts  further  doubt  on  the 
authenticity  of  many  of  these  models. 

VB.  Moment  Equation  Methods 

In  an  attempt  to  make  these  turbulent  energy  and  scale  models  more  reliable  and 
rigorous,  mathematical  approaches  have  attempted  to  generalize  the  intuitive  models  described 
above  into  more  formal  closure  procedures  (19,20,21,22).  The  physical  variables  are  decom¬ 
posed  into  a  mean  part  representing  the  macroscopic  flow  and  a  fluctuating  part  representing  the 
subgrid  scale  turbulence.  Equations  for  the  mean  and  fluctuating  components  are  derived  in 
terms  of  successively  higher  order  nonlinear  correlations  of  the  fluctuating  quantities.  The 
hierarchy  of  moment  equations  can  only  be  closed,  however,  by  making  an  assumption  that  the 
unknown  higher-order  terms  can  be  written  as  combinations  of  lower  order  quantities.  In  prac¬ 
tice,  !he  hierarchy  of  equations  is  closed  by  semi-empirical  arguments  which  range  from  very 
simple  guesses  for  an  exchange  coefficient  to  much  more  sophisticated  hierarchies,  in  which  the 
ultimate  closure  is  very  remote  from  any  physical  bases"  (8). 
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These  closure  approximations  often  contain  a  number  of  non-dimensional  parameters, 
which  are  hoped  to  be  universal  in  the  sense  of  being  constant  or  nearly  constant  over  a  wide 
range  of  turbulent  flows.  Of  course  this  approach  also  suffers  the  criticism  that  so  many  free 
parameters  become  potentially  available,  if  the  hierarchy  is  carried  far  enough,  that  any  func¬ 
tional  behavior  could  be  fitted.  Nevertheless,  for  certain  non-reactive,  incompressible, 
constant-density  flows,  satisfactory  agreement  between  the  models  and  selected  experiments  is 
obtained. 

The  extension  of  the  moment  methods  to  reactive  flows  has  used  the  idea  of  a  probability 
density  function  (PDF),  [23,24]  which  attempts  to  relate  time  averaged  macroscopic  quantities 
to  their  instantaneous  local  values.  Choice  of  the  PDF  is  close  to  arbitrary,  although  physical 
significance  is  being  attached  to  whether  delta  functions  or  trapezoids  are  used  [24].  A  major 
conceptual  problem  lies  in  the  fact  that  many  calculations  are  extremely  sensitive  to  the  form  of 
the  PDF. 

A  related  problem  centers  around  the  exact  order  and  form  of  closure.  It  was  originally 
hoped  that  by  going  to  higher  and  higher  order  models,  more  accuracy  would  be  obtained.  The 
current  feeling  is  that  the  simple  two-  or  three-equation  models  do  as  well  as  the  five,  six,  or 
twenty-six  moment  models.  Higher  orders  yield  more  free  parameters  which  can  be  fit  to  data. 
But  the  resulting  models  are  no  more  universal. 

Another  point  which  should  be  mentioned  here  is  basically  a  philosophical  problem.  The 
turbulence  modelling  schemes  currently  used  convert  the  detailed  consequences  of  microscopic 
motions,  which  are  thought  of  as  stochastic  random  phenomena,  into  absolutely  deterministic 
effects  driving  flow  on  the  macroscopic  scale.  Thus  in  a  reacting  fluid  system  which  is  modelled 
by  a  set  of  laminar  fluid  equations  and  a  prescription  for  the  effects  of  turbulence,  a  given  ini¬ 
tial  configuration  will  always  lead  to  the  same  final  state.  Even  though  the  real  fluid  system  is 


22 


NRL  MEMORANDUM  REPORT  4137 


being  described  in  the  model  by  only  a  small  finite  number  of  degrees  of  freedom,  the  answers 
predicted  by  the  model  are  completely  deterministic.  There  are  an  infinite  number  of 
configurations  of  the  real  system  which  all  would  be  represented  in  the  model  by  the  same  final 
representation,  and  each  of  these  real  configurations  should  really  lead  to  a  different  final  state. 
The  variance  of  these  final  states  within  a  given  finite  model  representation  has  to  be 
determined,  and  is  the  subject  of  extensive  mathematical  and  physical  research.  Distinct  classes 
of  solutions  seem  to  exist  in  simplified  physical  problems.  These  solutions  depend  to  whether 
the  problem  is  ergodic,  i.e.,  gives  rise  to  periodic,  quase-periodic,  or  non-periodic  solutions  [3]. 
Whether  this  mathematical  classification  bears  corresponding  physical  significance  is  not  deter¬ 
mined. 

It  is  by  no  means  clear  that  increased  study  and  development  of  these  mathematical 
phenomenological  approaches  is  the  way  to  go.  In  addition  to  the  problems  described  above, 
there  are  other  difficulties  as  well.  The  first  of  these  involves  modeling  turbulent  transients 
(bursts)  and  turbulence  onset.  In  most  physical  systems,  turbulence  arises  naturally  as  a  conse¬ 
quence  of  an  unstable  macroscopic  flow  configuration  and  dies  out  when  the  destabilizing  situa¬ 
tion  abates.  As  mentioned  in  Section  III,  specific  local  evaluation  of  stability  criteria  should  be 
a  part  of  any  reliable  model,  although  it  does  not  seem  to  be  included  in  any  of  the  moment 
models.  In  fact  the  description  of  turbulent  shear  flows  in  terms  of  a  mean  profile  in  at  least 
some  flows  (mixing  layer  and  boundary  layer)  has  been  shown  to  be  the  result  of  a  superposi¬ 
tion  of  nearly  deterministic  large-scale  structures  having  random  phases  [8].  The  classic  experi¬ 
ments  by  Roshko  and  coworkers  demonstrate  this  fact  clearly  (12,13). 

Another  aspect  of  these  transient-vs-steady-state  problems  is  shown  in  Figure  11.  Typical 
flow  streamlines  in  both  the  steady  state  and  the  transient  laminar  computations  of  flow  in  a 
simple  two-dimensional  combustor  are  illustrated.  The  lower  panel  shows  realistic  macroscopic 
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fluctuating  eddies  carrying  pockets  of  relatively  pure  oxidizer  to  the  top  of  the  chamber  and 
relatively  pure  fuel  to  the  bottom.  Any  phenomenological  turbulence  models  based  on  the 
steady  state  flow  of  the  top  panel  clearly  cannot  convect  pure  fuel  or  oxidizer  into  the  far 
corners  of  the  chamber,  for  it  would  involve  crossing  streamlines.  Diffusive  approximations 
will  also  not  be  able  to  get  the  pure  fuel  or  oxidizer  there  either.  Therefore,  suitably  modified 
approaches  to  the  modelling  are  required  to  correctly  treat  this  problem. 

Another  technical  issue  arises  from  the  overall  complexity  of  the  reactive  flow  problem. 
Current  turbulence  models  generally  include  approximations  developed  to  describe  incompressi¬ 
ble,  non-reacting  systems.  Energy  released  as  a  result  of  the  turbulent  mixing  and  the  subse¬ 
quent  chemical  reactions  which  occur  is,  as  described  extensively  in  Section  III,  qualitatively  a 
very  different  problem. 

VC.  A  Physically  Motivated  Phenomenological  Approach 

An  alternative  approach  to  the  multi-equation  and  multi-moment  approaches  involves 
actually  modelling  the  probability  distribution  function  (PDF)  based  on*  certain  physical 
processes  we  know  occur.  The  models  must  also  provide  procedures  for  deriving  the  time- 
averaged  fluid  quantities.  Methods  of  this  type  have  been  developed,  and  tested  by  a  number 
of  authors  (25,26,27,281.  Their  hope  is  that  by  building  some  of  the  physics  into  the  PDF,  the 
computed  mean  quantities  will  be  less  sensitive  to  those  of  its  feature  which  are  not  known. 

Spalding's  ESCIMO  theory  [28]  is  one  example  of  this  type  of  model  which  is  currently 
being  developed  Mixtures  of  different  composition  for  interleaved  folds  which  are  born, 
stretch  and  die  as  a  function  of  time.  Figure  12  shows  a  fold,  in  which  the  two  species  A  and 
B,  initially  have  sharp  gradients  across  their  interface.  The  gradients  are  smoothed  and 
broadened  by  molecular  diffusion.  Only  the  direction  normal  to  the  interface  is  considered,  so 


24 


NKI  MEMORANDUM  RII’ORI  4137 


multidimensional  effects  appear  as  compression  and  dilation  of  the  gradients  in  the  one  spatial 
dimension  which  is  represented.  The  internal  gradients  and  structures  of  the  folds  are  only 
driven  by  the  macroscopic  flow  but  not  derived  from  it  or  coupled  directly  to  it.  In  terms  of 
the  reactive  surface  concepts  introduced  in  Section  II  to  describe  turbulent  mixing  with  reac¬ 
tions,  the  ESCIMO  theory  concentrates  on  the  perpendicular  diffusive  mixing  and  associated 
detailed  kinetics  where  the  reactive  flow  is  laminar  on  the  small  scale.  The  effect  of  mutual 
interference  in  adjacent  folds  of  the  surface  are  neglected. 

This  theory  deals  with  the  evolution  of  a  fold  after  it  develops.  A  number  of  PDF’s  are 
defined  for  a  fold  which  are  functions  of  its  age,  thickness,  reactedness,  etc.  These  PDF’s  may 
then,  hopefully  be  used  to  obtain  averaged  fluid  quantities.  However,  one  encounters  a  funda¬ 
mental  problem  here:  there  is  no  way  of  knowing  what  the  size  distribution  is  for  these  folds  at 
their  birth.  What  we  have  is  a  number  of  small-scale  one-dimensional  laminar  problems,  which 
must  be  related  to  the  macroscopic  flow.  This  has  been  termed  the  "demography"  problem. 
Spalding  has  discussed  launching  a  large  number  of  these  independent  Lagrangian  reactor  par¬ 
cels  in  a  given  steady-state  flow  and  then  performing  statistical  analyses  on  the  ensemble  of 
these  parcels  from  their  "birth"  to  their  "death.”  Since  each  parcel  launched  this  way  is  an 
independent,  well  resolved  one-dimensional  reactive  flow  problem,  the  computational  cost 
would  be  enormous  even  in  simple  flows.  The  practical  cost  will  be  even  higher  because  the 
laminar  large-scale  flow  must  respond  to  the  short  scale  fluctuations.  Therefore,  the  "demogra¬ 
phy",  thought  of  as  an  a  posterior  averaging  procedure  should  be  performed  in  concert  with  the 
evolving  hydrodynamics  calculation  to  model  many  real  applications.  Again,  we  are  faced  with 
a  situation  where  we  hope  that  the  end  result  of  the  calculations  will  not  be  too  sensitive  to  the 
initial  choice  for  the  fold  distribution. 
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We  see,  therefore,  that  the  ESCIMO  approach  is  both  a  simplified  ab  initio  model  which 
can  teach  us  something  about  the  local  mixing  and  chemical  kinetic  aspects  of  turbulence  and 
may  be  the  basis  of  a  potentially  useful  phenomenological  model  of  turbulence.  Studies  of  the 
time-history  of  a  fold  may  include  effects  such  as  complex  chemical  reactions,  molecular 
diffusion,  etc.  It  is  one  way  to  study  localized  mixing  in  the  microscale  under  idealized  cir¬ 
cumstances.  The  model  is  limited,  however,  in  its  predictive  capability  by  those  physical 
processes  which  have  not  been  included,  as  well  a  by  the  open-ended  question  of  how  the  par¬ 
cels  interact  with  the  macroscale.  A  parameterized,  coupled,  dynamic  model  of  the  subcell  tur¬ 
bulence  is  still  required  for  use  in  detailed  modelling  calculations  which  is  capable  of  represent¬ 
ing  and  predicting  multiscale  phenomena  such  as  occur  in  instability  cascade. 

VD.  Local  Spectral  Dynamics 

Suppose  every  computational  fluid  element  has  local  subscale  motions  which  must  be 
described  simply  and  accurately.  Further,  suppose  these  local  motions  can  be  described  in 
important  respects  by  a  local  spectrum  of  fluctuations.  Then  an  efficient  and  flexible  turbulence 
model  might  be  developed  by  modelling  the  local  spectrum, 

Ek  =  \\PV?\. 

in  each  fluid  element  using  just  a  few  degrees  of  freedom.  Here  p  is  the  mass  density  and  Vk  is 
the  A:lh  wavelength  component  of  the  fourier  transform  of  the  fluid  velocity.  We  further 
assume  for  now  that  only  a  scalar  spectrum  is  needed  to  model  the  local  motions  and  that  the 
region  which  has  to  be  modelled  is  limited  as  shown  in  Figure  2.  For  wavenumoers  larger  than 
kdlsl,  viscous  dissipation,  thermal  conduction,  and  molecular  diffusion  effects  dominate  and  the 
fluctuation  spectra  should  be  predictable  given  the  values  at  kdlss.  For  wavenumbers  smaller 
than  kcfih  the  macroscopic  hydrodynamics  equations  resolve  these  scales  explicitly. 
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In  between  kdlSi  and  some  form  of  discrete  representation  of  the  localized  turbulence 
short  wavelength  spectrum  is  needed.  The  best  form  has  not  yet  been  determined.  In  the 
figure  three  discrete  values  of  the  spectrum  at  k<,  k2 ,  and  k2  are  shown  breaking  the  turbulent 
interactions  regime  into  three  finite  width  cells  or  bins.  The  optimal  number  of  such  cells  must 
also  be  determined.  At  least  three  bins  are  necessary  to  allow  for  transient  instability  onset, 
changing  turbulent  length  scales,  and  multiscale  phenomena.  More  bins  will  undoubtedly 
improve  flexibility  and  accuracy  of  the  model.  In  earlier  idealized  tests  of  this  approach  on  the 
Kelvin-Helmholtz  instability,  it  was  necessary  to  use  at  least  two  bins  per  decade  of  characteris¬ 
tic  scale  length. 

To  advance  Ek ,  Eki,  Ek} ,  the  localized  spectral  strengths  at  k\,  k2,  k 3,  a  set  of  ordinary 
differential  equations  is  written  down  which  are  somewhat  analogous  to  the  coupled  rate  equa¬ 
tions  for  chemical  kinetics.  For  each  macroscopic  cell,  there  are  a  series  of  equations  for  the 
form 

dEk 

~~r-  “  £  «?*,<>;)  "  Lt(yj)) 

ai  j 

where  Qk  and  Lk  are  production  and  loss  terms  for  spectral  cell  k,.  The  (y,)  refer  to  the  par¬ 
ticular  processes  which  can  act  as  a  production  ( Q )  or  loss  (L)  terms  and  they  represent  the 
hydrodynamic  channels  for  turbulent  energy  reapportionment,  chemical  reactions,  diffusion 
processes,  etc.  These  source  and  sink  terms  couple  the  various  wavelength  bins  in  much  the 
same  manner  as  specific  chemical  reactions  or  reaction  types  couple  different  chemical  kinetic 
species.  All  the  conservation  conditions  on  total  momentum  and  energy  in  the  spectrum  apply 
with  the  added  complication  of  source  and  sink  terms  at  both  long  and  short  wavelength. 
Energy  dissipated  at  short  wavelengths  appears  as  localized  heat  in  the  macroscopic  energy 
equation  and  as  turbulent  viscosity  in  the  macroscopic  momentum  equation. 
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In  preliminary  studies,  the  obvious  anisotropy  of  some  of  the  hydrodynamic  channels  will 
be  ignored,  and  it  is  hoped  that  the  use  of  several  spectra  will  not  be  necessary.  Fully  cali¬ 
brated  and  theoretically  understood  parameterizations  of  the  reactive  surface  area  evolution  as  a 
function  of  spectral  structure  (Section  II)  will  be  needed  if  accurate  chemical  kinetics  models 
are  to  be  included.  Further,  the  complementary  results  of  a  complete  set  of  ESCIMO-like 
Lagrangian  reactive  flow  calculations  will  need  to  be  folded  in.  Certainly  only  a  careful,  quanti¬ 
tative  analysis  of  the  specific  hydrodynamic  channels  and  their  interactions  can  point  the  way  to 
valid  simplifications  in  this  miserable  modelling  morass.  This  analysis  must  provide  a  way  to 
represent  the  basic  fluid  instabilities  as  source  and  sink  terms,  probably  in  terms  of  instability 
growth  rates  derived  from  separately  performed  ab  initio  or  microscopic  calculations.  Consider¬ 
able  success  with  a  variant  of  this  approach  has  been  achieved  in  a  number  of  plasma  physics 
problems  where  small-scale  micro-instabilities  affect  macroscopic  transport  [9,10,29,30,31],  A 
similar  approach  has  been  postulated  for  nonreactive  flow  by  Brodkey  [32].  However,  his 
source  and  sink  terms  did  not  explicitly  identify  individual  hydrodynamic  channels  or  a  way  of 
incorporating  the  effective  boundary  conditions  of  the  macroscopic  gradients  on  the  subscale 
quantities.  Studies  of  the  effect  of  chemical  energy  release  on  the  spectrum,  £*,  have  been 
made  by  Eschenroeder  [33]. 

Currently  this  local  spectral  dynamics  approach  is  being  pursued  in  the  Laboratory  for 
Computational  Physics  at  NRL.  The  goal  is  to  use  a  combination  of  ab  initio  calculations  and 
analytical  theory  to  provide  insight  into  the  form  of  the  Pk  and  Lk  functions.  This  approach  is 
appealling  physically  because  it  allows  all  the  important  aspects  of  the  six  criteria  in  Section  III 
to  be  met.  The  model  may  include  the  growth  and  decay  of  the  important  hydrodynamic  chan¬ 
nels  and  thus  potentially  allows  for  onset  and  transient  phenomena.  Chemical  reactions  may  be 
included,  at  the  very  least  through  an  enhancement  of  the  spectrum  at  wavelengths,  between 
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kctu  and  kdm.  Since  we  clearly  must  accept  a  prediction  which  is  true  only  on  average,  the 
neglect  of  all  but  a  few  of  the  internal  turbulent  fluid  degrees  of  freedom  in  each  cell  ensures 
this  and  provides  efficiency.  Scaling  is  ensured  through  a  source  term  imposed  on  the  equa¬ 
tions  at  kce„.  Thus  the  procedure  sketched  out  above  does,  at  least  in  principle,  meet  the  six 
conditions  postulated  earlier  as  a  necessary  if  not  sufficient  for  a  reactive  flow  turbulence  sub¬ 
model. 
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VI.  Conclusion 

Computational  simulation  techniques  play  two  roles  in  reactive  flow  turbulence  and  tur¬ 
bulent  mixing  research.  Ab  initio  calculations  of  basic  turbulence  processes,  in  which  the  fluid 
dynamic  and  chemical  phenomena  are  coupled  on  the  microscale  must  be  performed  to  extend 
our  understanding.  In  addition,  macroscopic  phenomenological  models  of  turbulent  mixing  are 
required  which  are  efficient  enough  for  regular  use  in  practical  applications.  These  two  roles 
complement  each  other  and  hence  must  be  closely  coupled  in  a  comprehensive  approach.  Only 
the  ab  initio  calculations  can  bring  out  new  knowledge  about  turbulence  itself.  The 
phenomenological  models  contain  no  more  physics  than  is  built  into  them  and  thus  are  wholly 
dependent  upon  theory,  experiment,  and  ab  initio  calculations  for  their  content. 

The  obstacles  confronting  the  detailed  modeller  who  must  include  the  effects  of  turbulent 
mixing  depend  on  whether  an  ab  initio  or  a  phenomenological  approach  is  adopted.  Using  ab 
initio  calculations,  the  modeller  deals  with  a  straight-forward  representation  of  the  fluid  dynam¬ 
ics,  but  adequate  resolution  of  varied  space  scales  will  always  be  a  major  stumbling  block. 
Improved  computers  and  numerical  algorithms  help,  but  for  years  to  come,  detailed  simulations 
will  be  confined  to  idealized  problems  in  simplified  geometries.  The  mtuor  hurdle  to 
phenomenological  modelling  is  the  choice  and  implementation  of  a  suitable  representation.  A 
successful  separation  and  then  asymptotic  coupling  of  the  macroscopic  and  microscopic  scales 
makes  the  resolution  problem  becomes  less  important. 

In  this  paper,  the  physical  origins  of  turbulence  and  the  desirable  properties  of  a  general 
time-dependent  subgrid  turbulence  model  have  been  described.  The  localized  and  transient 
aspects  of  turbulent  mixing  have  been  stressed  along  with  the  dynamics  of  reactive  interfaces 
and  the  strong  effects  of  chemical  energy  release.  It  is  hoped  that  by  thoroughly  understanding 
these  aspects  of  turbulent  mixing  in  reactive  flows  a  suitable  phenomenological  model  can  be 
built  which  displays  the  desired  properties. 
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F,gure  I  -  A  schematic  diagram  showing  the  in.eract.ons  between  the  three  different  aspects  of  reactive  flow  mixmg 
Lentdied  as  1  .hydrodynamic  channels  excited  by  velocity,  density,  and  pressure  grad.ents.  2  react, ve  m.e  face 
dynamics  where  molecular  interpenetration  of  reactive  species  is  enhanced,  and  3)  detailed  chemical  kinetics  who.  e 
pansion  and  other  hea,  release  effects  feed  energy  back  into  the  fluid  dynamics  to  drive  additional  hydrodynamic  chan- 


k  <  kcen.  the  transient,  macroscopic  fluid  calculation  may  be  performed  For  short  wavelengths,  k  >  diffusion 
and  dissipation  dominate  so  simplified  analyses  suffice  The  intermediate  region,  called  inertial,  subgrid,  or  turbulent, 
is  broken  into  bins  centered  at  discrete  wavenumbers  k;.  k ,  The  hope  is  that  an  efficient,  tractable  set  of  equa¬ 
tions  can  be  written  down  for  the  spectral  strength  of  the  turbulence  at  each  wavenumber  which  matches  suitably  to 
the  two  known  regions  to  the  right  and  left 
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Figure  3  —  The  active  role  of  expansion  in  the  turbulence  of  reactive  flows.  A  sharp  density  gradient  is  shown  in  the 
vicinity  of  a  molecularly  mixed  pocket  of  fuel  and  oxidizer  where  reactions  and  heat  release  are  occurring.  The  expan¬ 
sion  induced  pressure  gradients  interact  with  the  density  gradient  to  produce  vorlicity  in  the  flow. 
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Figure  4  —  The  passive  role  of  expansion  in  the  turbulence  of  reactive  flows.  Here  previously  mixed  pockets  of  fuel 
and  oxidizer  have  reacted  and  thus  expanded  to  low  density.  These  density  gradients  again  generate  vorticity  on  a  local 
small  scale  as  they  interact  with  global  pressure  gradients  resulting  from  larger  scale  circulations  of  the  fluid. 
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(upper  left),  the  surface  separating  A  and  B  is  shown  as  flat  and  two-dimensional  in  a  three-dimensional  volume  of 
fluid.  At  later  times  ,  (upper  right)  and  2  (lower  left),  the  surface  progressively  deforms  as  the  result  of  the  local 


Lagrangian  motion  of  surface  points  shown  as  small  arrows  to  indicate  local  direction  of  flow. 
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Figure  7  —  Two  different  cases  for  reactive  interface  area  increase  are  shown.  In  the  upper  right  a  relatively  flat  fluffy 
surface  results  from  a  turbulent  velocity  spectrum  enhanced  at  short  wavelengths  where  molecular  diffusion  is  strong. 
Although  the  ideal  surface  area  should  increase  exponentially  as  shown  in  the  lower  right,  the  effective  surface  area  will 
increase  more  slowly  as  indicated  by  the  —  —  line  because  diffusion  can  nullify  small  wavelength  area  increases. 
When  diffusion  is  small  and  the  velocity  spectrum  large  at  long  wavelengths,  smooth  bulges  result  as  shown  in  the 
lower  left  Exponential  area  growth  is  possible  but  at  a  lower  rate  in  this  case. 
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t  =  0.0  sec  t  =  0.3  sec 


density  contour  evolution  of  initial  checkerboard 


Figure  8  —  Contour  plots  of  density  during  a  detailed  two-dimensional  reactive  flow  simulation  designed  to  show  the 
effects  of  density  differences  on  the  generation  of  short  wavelength  turbulence  and  mixing.  Four  regions  of  alternating 
high  and  low  mass  density  are  established  in  a  checkerboard  pattern  at  lime  /  -  0.0  sec  (see  upper  left  panel  of  Fig  9). 

Contours  at  p„  +  «(p6  -  p„)  are  shown  where  <r  -  0.1,  0,  2 . 0.9.  After  0.3  seconds  the  beginning  of  convective 

mixing  is  seen  due  to  an  initially  impressed  random  velocity  field.  At  the  two  later  times  shown  (lower  panels),  molec¬ 
ular  mixing,  enhanced  by  short  wavelength  convection,  leaves  only  the  contours  at  «  -  0.4.  05,  and  0.5. 
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Ab  Initio  Simulation  of  Reactive  Flow  Mixing 


t  =  0.0  sec 


t  =  0.3  sec 


density  surface  evolution  of  initial  checkerboard 

f  igure  •>  -  Plots  of  density  us  .1  function  of  v  ami  1  al  four  times  for  the  ah  intro  reactive  flow  nosing  calculation  of 
Fig  X  As  can  he  seen,  molecular  mixing  rapidly  smooths  out  the  jagged  interfaces  which  result  from  the  short 
wavelength  components  of  the  "turbulence"  field  at  the  beginning  of  the  calculations 
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and  developing  phenomenological  subgrid  models. 


THE  PROBLEM  WITH  STEADY  FLOW  APPROXIMATIONS 
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Figure  1 1  —  A  practical  and  philosophical  problem  with  slcady  flow  approximations  The  upper  panel  shows  an  ideal¬ 
ized  flow  assumed  to  be  the  steady  state  time  average  of  a  turbulent  mixing  flow  Pure  oxidizer  and  pure  fuel  cannot 
cross  the  separating  streamline  so  pure  oxidizer  in  the  upper  right  corner  of  the  combustor  or  pure  fuel  in  the  lower 
left  corner  is  an  impossibility.  In  the  real  fluctuating  flow  of  the  lower  panel,  whorls  and  vertices  ensure  that  there  are 
large  scale  variations  in  flow  and  that  fuel  and  oxidizer  are  carried  intermittantly  to  the  far  corners  of  the  system 


Figure  12  —  A  global  perspective  on  the  ESCIMO  approach.  Each  element  of  the  evolving  reactive  surface  between 
two  species  can  be  viewed  as  supporting  a  one-dimensional  reactive  flow  problem  where  only  the  spatial  direction  nor¬ 
mal  to  the  surface  is  considered.  Rather  than  emphasizing  the  global  dynamics  of  the  reactive  surface,  this  approach 
focuses  on  the  complementary  reaction  kinetics  of  the  diffusively  interpenetrating  fluids  across  the  interface  as  shown 
in  the  lower  right  Both  approaches  are  crucial  to  extending  our  understanding  of  reactive  (low  mixing 


